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Nonlinear periodic systems, such as photonic crystals and Bose-Einstein condensates (BECs) 
loaded into optical lattices, are often described by the nonlinear Schrodinger/Gross-Pitaevskii equa- 
tion with a sinusoidal potential. Here, we consider a model based on such a periodic potential, with 
the nonlinearity (attractive or repulsive) concentrated either at a single point or at a symmetric set 
of two points, which are represented, respectively, by a single 5-function or a combination of two 
^-functions. With the attractive or repulsive sign of the nonlinearity, this model gives rise to ordi- 
nary solitons or gap solitons (GSs), which reside, respectively, in the semi-infinite or finite gaps of 
the system's linear spectrum, being pinned to the <5-functions. Physical realizations of these systems 
are possible in optics and BEC, using diverse variants of the nonlinearity management. First, we 
demonstrate that the single 5-function multiplying the nonlinear term supports families of stable 
regular solitons in the self-attractive case, while a family of solitons supported by the attractive 
(5-function in the absence of the periodic potential is completely unstable. In addition, we show 
that the 5-function can support stable GSs in the first finite bandgap in both the self-attractive and 
repulsive models. The stability analysis for the GSs in the second finite bandgap is reported too, 
for both signs of the nonlinearity. Alongside the numerical analysis, analytical approximations are 
developed for the solitons in the semi-infinite and first two finite gaps, with the single 5-function 
positioned at a minimum or maximum of the periodic potential. In the model with the symmetric 
set of two 5-functions, we study the effect of the spontaneous symmetry breaking of the pinned soli- 
tons. Two configurations are considered, with the 5-functions set symmetrically with respect to the 
minimum or maximum of the underlying potential. 
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I. INTRODUCTION 

In the course of the last decade, a great deal of interest has been drawn to theoretical and experimental studies 
of the nonlinear dynamics in systems with periodic potentials. Physical realizations of this topic are well known in 
nonlinear optics [TJ [3] and Bose-Einstein condensates (BECs) . In optical media, the periodic (lattice) potentials may 
be created as permanent or virtual ones (in the latter case, these are photonic lattices induced in photorefractive 
crystals [2111]). In BECs, similar potentials can be induced in the form of optical lattices (OLs), i.e., interference 
patterns formed by laser beams shone through the condensate. The OLs make it possible to study a great variety of 
dynamical effects in BECs [21 [5] . Similar periodic potentials may also be imposed by magnetic lattices [5] . 

It is commonly known that, in the free space, stable one-dimensional (ID) solitons exist in optical waveguides and 
BECs with the attractive cubic nonlinearity, while in the 2D and 3D geometry the solitons are unstable to the collapse 
[7j. On the other hand, it has been predicted that 2D [SHU] and 3D [PJH2] solitons can be stabilized by dint of the 
corresponding OLs. Moreover, low-dimensional OLs, i.e., quasi-lD and quasi-2D lattices in the 2D [T3] and 3D (T31[T3] 
space, respectively, also provide for the stabilization of fully localized multidimensional solitons. A related prediction 
is the existence of stable 2D [14] [15] and 3D [16] solitons in models with radial OLs. 

In BECs with repulsive interactions between atoms, i.e., the repulsive intrinsic nonlinearity, solitons cannot exist 
in free space, but gap solitons (GSs) may be supported by the OL. The principle behind the formation of the GS is 
that the periodic potential can invert the sign of the effective mass of collective excitations, which may then balance 
the repulsive nonlinearity. In the ID case, several species of GSs are known, which include fundamental solitons and 
their two- and three-peak bound complexes, in the first and second bandgaps [XT] [15], and subfundamental solitons 
in the second gap (the latter means a twisted soliton squeezed into a single cell of the potential lattice, whose norm 
is smaller than the norm of the fundamental soliton existing at the same value of the chemical potential) [T9lI21j . 
The origin of the GS families may be traced back to bifurcations generating them from Bloch waves at edges of the 
bandgaps (151123]. 

Multidimensional GSs are represented by fundamental solitons [3^1125] and gap-type vortices [321 [35] I39TI3T] . In 
particular, the simplest gap vortices are composed of four density peaks, and fall into two different categories: densely 
packed squares (alias off-site-centered vortices), in which the center is positioned around a local maximum of the OL 
potential [24, 29, 31 , and rhombic (on-site-centered) configurations, featuring a nearly empty lattice cell in the middle 
[33H36] . While GS families in 2D are also generated by bifurcations from the respective Bloch waves [23] [27] , there 
are gap- vortex families which do not originate from such bifurcations [28] , 

In the experiment, an effectively one-dimensional GS, composed of few hundred atoms, was created in the condensate 
of 87 Rb [32]. In another experiment, which involved a stronger OL, broad confined states were created [33]. They 
were interpreted as modes intermediate between the GSs and extended nonlinear Bloch waves |34j . 

In a different framework, solitons can be supported by a spatial modulation of the local nonlinearity (a particular 
case of the "nonlinearity management" |35j). In terms of solid-state physics, an effective potential structure induced 
by means of this method is called a pseudopotential [36] . Realizations of such structures are possible in optics and 
BECs, where they are expected to give rise to many new effects, see recent review [37 . In optical media, nonlinearity- 
modulation profiles can be induced by means of nonuniform distributions of resonant dopants, similar to those in 
resonantly-absorbing Bragg reflectors, see original works [35] and review [3S]. Such dopant-density patterns can be 
created by means of available technology [10]. For matter waves in BECs, similar nonlinear profiles can be created via 
the spatial modulation of a s (x), the local value of the s-wave scattering length. Such modulations can be induced by 
means of the Feshbach resonance, controlled by a nonuniform dc magnetic field [21] or resonant optical field [321 143j . 
It was also predicted that the Feshbach resonance may be controlled by dc electric fields [22] ■ Note that the Feshbach 
resonance allows one to create a pseudopotential corresponding to a sign-changing function a s (x) (actually, this is 
possible in the condensate of 7 Li [45]). which implies the spatial alternation between the attractive and repulsive signs 
of the nonlinearity. 

In the context of the BECs, many works dealt with solitons and related dynamical states in the framework of ID 
pseudopotential structures [46] . Nonlinear structures based on the spatial modulation of the Kerr coefficient were also 
considered in optics, giving rise to similar states [27] • Solitons supported by 2D pseudopotentials were theoretically 
studied too, with a conclusion that it is much more difficult to stabilize them than using linear OL potentials [48] . 
The results accumulated in the studies of solitons in nonlinear lattices are summarized in a recent review [37] . 

A specific example of such nonlinear pseudopotential settings is a model where the nonlinearity is concentrated 
at a single point, which is represented by a (^-function. A prototypical model of this type was introduced in Ref. 
[49j . It may represent a planar linear waveguide with a narrow nonlinear stripe embedded into it (for the case of the 
second-harmonic-generating quadratic nonlinearity, spatial solitons in a similar setting were considered in Ref. (50) ) . 
In terms of the BECs, localized nonlinearity may be induced through the Feshbach resonance imposed by a focused 
laser beam. It should be said that a more realistic model would also include a uniform background nonlinearity, as it 
is difficult to create a setting where it might be completely eliminated. Nevertheless, it makes sense to consider, as 
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the basic one, the model in which the nonlinearity is represented solely by the localized terms, as it is reasonable to 
assume that the weakly nonlinear background will not lead to a drastic change of the results. Extending the studies 
in this direction, a model of a double-well pseudopotential, based on a symmetric set of two ^-functions, or their 
regularized counterparts, was introduced, with the purpose of studying the spontaneous symmetry breaking (SSB) 
of localized modes [5T ]. In particular, it is possible to find full analytical solutions for symmetric, asymmetric and 
antisymmetric states in the model with two ideal 5-functions. In this model, where the SSB bifurcation is of the 
subcritical type, the symmetric solutions are stable up to the bifurcation point. Beyond this point, the symmetric 
states and the emerging asymmetric states are unstable, as well as all the antisymmetric ones. Symmetry breaking 
in a circular nonlinear lattice, with a smooth spatial modulation of a s , was studied, in the framework of both the 
Gross-Pitaevskii equation (GPE) and many-body quantum system, in Ref. |52j . 

A natural extension of the settings outlined above is a model integrating the linear OL potential and the nonlinearity 
concentrated at one or two points, represented by the respective 5-functions or their regularized versions. These 
systems are the subject of the present work. In particular, we demonstrate that, while cusp-shaped solitons pinned 
to the 5-function multiplying the attractive cubic nonlinearity turn out to be unstable, the linear periodic potential 
readily stabilizes them. Another issue of obvious interest is whether the strongly localized repulsive nonlinearity may 
support GSs in finite bandgaps of the OL- induced spectrum (we demonstrate that this is possible indeed). 

The rest of the paper is organized as follows. The model is formulated in section [TTJ For the system including 
the single 5-function and OL potential, analytical approximations, based on the perturbation theory, are developed 
for the pinned modes in the first and second finite bandgaps, as well as in the semi-infinite gap, in section |III| A 
comparison with numerical findings is performed too. Detailed results of the numerical analysis for the existence and, 



most important, stability of the pinned modes in the semi-infinite, first and second gaps are reported in section IV 



Both the attractive and repulsive signs of the 5-functional nonlinearity are considered, for different positions of the 
5-function with respect to the underlying lattice. In section [V| the analysis is reported for symmetric, antisymmetric 
and asymmetric modes supported by a pair of the 5-functions, positioned symmetrically with respect to a maximum 



or minimum of the OL potential. The paper is concluded by section VI 



II. THE MODEL 



The model featuring the cubic nonlinearity represented by the 5-function was introduced in Ref. 49J , in the context 
of tunneling of interacting particles through a junction: 

ii>t + \^xx ~ crS(x)\ip\ 2 ip = 0, (1) 

where ip is the mean- field wave function in the BEC, or the local amplitude of the guided electromagnetic field in the 
context of optics (in the latter case, time t is replaced by the propagation distance), with a = +1 and —1 corresponding 
to the repulsive and attractive nonlinearity, respectively. Obviously, Eq. ([I]) amounts to the simple linear equation 
valid at x < and x > 0, iip t + [\/2')^j) lcx = 0, which is supplemented by the derivative-jump condition at x = 0, 
produced by the integration of Eq. (JlJ over an infinitely small vicinity of x = 0: 

tP x (x = +0)-i> x (x= -0) = 2a\4,(x = 0)\ 2 i>(x = 0). (2) 

Stationary states are looked for in the ordinary form, ip(x, t) — e~ v * <fi(x) , where \i is the chemical potential in BEC 
(— yU is the propagation constant in optics), and real function ip obeys equation 

^ + (1/2)0" - <t6(x)c/) 3 = 0. (3) 

An exact soliton solution to Eq. (|3| with a = — 1 is obvious: 

0o (x) = AeS=™*\, A 2 = v/^V (4) 

The norm of this solution, 

/+oo 
Mx)\*dx, (5) 
-oo 

is independent of fi, taking a constant value, N = 1, hence the formal application of the Vakhitov-Kolokolov (VK) 
criterion, dN/d^i < [71 153j. predicts neutral stability. In fact, numerical simulations of the evolution of pinned 
solitons Q with small perturbations added to them demonstrate a strong instability (not displayed here in detail): 
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the soliton either decays or suffers the collapse (formation of a singularity), if the exact norm of the perturbed soliton 
is, respectively, N < 1 or N > 1. These features, including the degeneracy of the norm and the instability leading 
to the decay or collapse, resemble those known for the 2D Townes solitons [7] in the free 2D space with the cubic 
attractive nonlinearity, or in the ID space with the quintic nonlincarity [5]. Moreover, it is easy to find a family of 
exact analytical solutions to Eq. (fib , with a = — 1, which explicitly describe the approach to the collapse at t — > CP: 



ip{x,t) 



Xq 



exp 



IX ) 



2t 



(6) 



where xq is an arbitrary real positive constant, the solution being valid at t < 0. Decaying solutions, at t > 0, are 
described by the same solution ([6]), with — xq replaced by xq > 0, the decay taking place at t — >• +oo. The norm of 
solution ^ is exactly N = 1 (irrespective of the value of xq), i.e., the same as that of stationary solution Q. 

In this work we introduce a natural extension of Eq. adding to it the periodic OL potential, with the objective 
to stabilize the solitons: 



ifpt + \^xx + £COs(2x)ip - aS(x - £)\i/}\ 2 i/> = 0. 



(7) 



Here, the OL potential, whose period is normalized to be tt, is V(x) — — ecos(2x). In its first period, < x < n, 
the minimum and maximum of the potential are located, respectively, at x — and x = n/2 (and vice versa for 
e < 0), while the 5-function is set at x = £, which does not necessarily coincide with the minimum or maximum of 
the potential. 

A modification of the model, which is considered in section [Vj deals with the nonlinearity represented by a pair of 
((-functions, placed symmetrically with respect to the minimum or maximum of the periodic potential: 



itpt + ~^xx + £ cos(2a;)V> - a [5(x - £) + 8(x + £)] |V>I V = 



(8) 



(recall the same double-delta nonlincarity, but without the OL potential, and solely with the attractive sign of the 
nonlinearity, a = — 1, was introduced in Ref. [ST]). In this case, the asymmetry measure of stationary modes is 
defined as 



f + °° \i,{ X t dx r: x Mx)? dx N+ N _ 



N 



(9) 



Dynamical invariants of Eqs. ^ and {[sj) arc norm N, given by Eq. and the Hamiltonian (written here for the 
latter equation), 



H = 



+ 00 " 



- \ip x \ 2 - ecos (2x) \i})(x)\ 2 



dx 



(10) 



+,- 



N is proportional to the number of atoms trapped in the BEC, or the total power of the trapped beam in optics. 
Control parameters of the models are e and £, along with N. 

As said above, in models with periodic potentials solitons may exist in bandgaps of the spectrum of the linearized 
version of the equation - either the semi-infinite gap or finite ones, separated by Bloch bands [4j [T7J [18] . For the 
models under consideration, the spectrum is the same as in the classical Mathieu equation (it is actually displayed by 
means of shaded areas in Fig. [l] below). 



III. PERTURBATION ANALYSIS OF THE SINGLE-DELTA MODEL 



Equation Q, as well as the corresponding linear Mathieu equation, can be treated by means of the perturbation 
theory if the OL strength, e, is a small parameter. Here, we consider the case when the ^-function is placed symmet- 
rically with respect to the periodic potential, i.e., at £ = 0, the objective being to construct approximate analytical 
solutions for GSs (gap solitons) supported by the repulsive [a = +1) or attractive (a = —1) point-wise nonlinearity, 
in the first and second finite bandgaps, as well as for ordinary solitons in the the semi-infinite gap. The respective 
form of the stationary equation is 

^ + (1/2)0" + [ecos(2a;) - o5{x)<j) 2 } = 0, (11) 



cf. Eq. (3). 
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A. The first finite bandgap 

1. Solutions of the linear equation 

In the case of small |e|, the first finite bandgap occupies a narrow interval of values of \i around /i = 1/2. Accordingly, 
wc set 



H=l/2 + u, with \v\ < 1. 
Then, approximate solutions corresponding to the GS in the first finite bandgap can be sought for as 



4>(x) 



[Acos(x) + B sin (\x\ 



(12) 



(13) 



where Ai is assumed to be a small positive coefficient (subscript 1 indicates that Ai pertains to the first finite bandgap). 

The substitution of ansatz (131 into Eq. (Ill, off point x — (which corresponds to the linear equation), and the 
subsequent analysis following the usual perturbation theory for the Mathieu equation (in other words, the asymptotic 
analysis of the linear parametric resonance [54] ) . yield the following equations for amplitudes A and B: 



A-A 1 S = ! 
B + X 1 A = 0. 



The solvability condition for the linear homogeneous system ( 14 1 is 

■ > 



In the first approximation, A^/2 may be neglected in the parentheses, which yields 



A a 



hence the solution exists for 



M < kl/2. 



(14) 



(15) 



(16) 



(17) 



In fact, Eq. (17) is the prediction provided by the perturbation theory for the width of the first finite bandgap, see 



Fig. 1(a) Further, it follows from Eqs. (14) that, in the first approximation 



(18) 



where condition (17) is taken into regard, to identify the correct sign. 



2. The nonlinear part of the solution 
For the stationary solutions, condition ([2| takes the form of 



A I ^ ) | x=0 = 2atp 1 , 1 1 



(19) 



where A stands for the jump of the derivative at x = 0. The substitution of Eqs. (13), (16), and (18) into Eq. (19) 
yields 



(20) 
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On the left-hand side of Eq. (20), the term \/s 2 /4 
\J e 2 /4 — v 2 is small, while + v) I (e/2 - v) 



2 is much smaller than the following term [in the general case, 
is not; in special cases, near the upper edge of the bandgap for 
•C |e/2 — v\, the term \J (e/2 + v) / (e/2 — v) becomes small, but 



e > 0, or the lower one for e < 0, with |e/2 H 
•\/e 2 /4 — ^ 2 is then still smaller]. Thus, Eq. (20) yields, for nontrivial solutions (A 2 ^ 0) 



A 2 ~ a sgn (e) 



/e/2- 



e/2- 



(21) 



This solution exists provided that sgn (e) = a, and it does not exist in the opposite case. In other words, it exists if 
the repulsive (5-function, with a = +1, is set at the local minimum of the OL potential, or the attractive (5-function is 
placed at the local maximum. Note that, although the result was obtained by means of the perturbation theory, the 
amplitude given by Eq. (21) is not small (which does not invalidate the pcrturbative treatment). 



In the first approximation, the calculation of the norm of the weakly localized soliton solution based on the above 
formulas yields 



(2Ai) _1 (A 2 + B 2 ) 



(22) 



In the case of the attraction (then, e is negative, as shown above), relation (22) yields 

3 < °' 



dN x _ dNi _ 
dn ~ dv (e/2-i'Y 



(23) 



hence, according to the VK criterion, the so predicted family of the GSs might be stable. Nevertheless, it is actually 
unstable, as shown below. As for the repulsive nonlinearity, the VK criterion is irrelevant for it (being sometimes 
replaced by an "anti-VK" condition [55]). 

The calculation of Hamiltonian (|10| in the present approximation yields 



H 



(24) 



[in this approximation, H is dominated by the gradient term in Eq. ( |10[ )]. The fact that this expression is always 
positive demonstrates that the GS cannot realize a ground state (it cannot correspond to an absolute minimum of 
the energy, which must be either negative or zero). Nevertheless, this does not mean that these solitons cannot be 
stable against small perturbations (actually, they may be metastable states, in comparison with the ground state). 



3. Comparison with numerical results 



The first finite bandgap, as predicted by the perturbation theory [Eq. (17)], is shown in Fig. l(a)| together with 
the numerically found borders of the first bandgap. In accordance with the above prediction, the soliton in the model 
with the repulsive or attractive nonlinearity exists only for e > and e < 0, respectively (this conclusion is true for 
all values of e and /j, in the first bandgap, not only for small e). 

As concerns the stability, the GSs supported by the repulsive (5-function placed at the minimum of the OL potential 
(er = +1, e > 0) are stable, while all the solitons generated in the finite bandgap by the attractive nonlinearity 
(cr = — 1, e < 0) are unstable [on the contrary to the formal prediction of the VK criterion, see Eq. (23)]. The 
detailed stability analysis is reported in the next section. 

Typical profiles of the soliton solutions are shown in Figs. [2] and [3j If e is not too large, the approximation quite 
accurately predicts the shape of the soliton. In particular, the analytical and numerical results agree very well for \x 
taken near the middle of the bandgap. Close to the upper edge (for e > 0) or lower edge (for e < 0), the analytical 
shape of the soliton is less accurate, as seen in Figs. |2(c)|(f ) and 3(a)|(d) In fact, this inaccuracy originates from the 
use of approximation (21 ) instead of the more general equation (20). 



Figure [4] shows the norm of the solitons versus /j,, for both negative and positive e, which correspond to the attractive 
and repulsive nonlinearities, respectively. As expected, the results become more accurate as |e| diminishes, while /i 
takes values farther from the upper or lower edges of the bandgap, for e > and e < 0, respectively. 
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FIG. 1. Borders of the first (a) and second (b) finite bandgaps (labeled as Gap I and Gap II, respectively), as predicted by 
the perturbation approximation, see Eqs. (17 1 and (32 l, in comparison with the numerically constructed bandgap structure, 



in which Bloch bands are shaded. The semi-infinite gap is represented by the white areas at the bottom of the panels. 



B. The second finite bandgap 

The perturbation theory for the Mathieu equation may also uncover the second finite bandgap, when 

H = 2 + u, with \u\ < 1, (25) 



cf. Eq. (12). To this end, an approximate solution to the linear equation is looked for as 

(j)(x) = e~ X2lxl [A + A 2 cos (Ax) + B 2 sin (A\x\) + A x cos (2x) + B x sin (2\x\ 



(26) 



neglecting terms ~ A2, we find, in the zeroth-order approximation, 

A, 



cf. ansatz (j 1 3|) adopted in the first finite bandgap. Substitution ansatz (|26|) into Eq. (Ill, off point x = 0, and 

(27) 



-Ay, A 2 = —Ay, B 2 = —By 

4 12 12 



Next, the homogeneous system of equations for Ay and By takes the following form, cf. Eqs. (14) derived in the first 
finite bandgap: 



vAy + -XjAy + e ( A) + ~A 2 } - 2\ 2 By = 0, 



uBy + -X 2 2 By + -eB 2 + 2X 2 Ay = 0, 



(28) 
(29) 



or, on substituting expressions (27), 



,-1 + 1)^-2^=0, 

2 \ 2 s 



(30) 



Like in the case of Eqs. (14), terms \ 2 /2 in the parentheses may be neglected in the lowest approximation, which 
yields 



(31) 
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(d) (e) (f) 



FIG. 2. Comparison between the analytical approximation (dashed lines) and the numerically found profiles (solid lines) for 
solitons in the first finite bandgap of the repulsive model, for e = 0.5 (a)-(c) and e = 1 (d)-(f). The results are shown for 
selected values of [i (indicated in each panel), close to the center of the gap or near each of its edges. 



As follows from this expression, the perturbation theory predicts the following form of the second finite bandgap, see 



Fig. 1(b) 



1/24 < v = v/e 2 < 5/24, 



(32) 



cf. Eq. (17) for the first finite bandgap. Then, in the lowest approximation, the relation between Bi and Ai is [cf. 
Eq. PIP 



(5 £ y 2 4)- v 

2A 2 



(33) 



Further, the jump condition at point x — keeps the form of Eq. ( 19 ), and, in the lowest approximation, the jump 
of the first derivative is dominated by term Bi sin (2\x\) in ansatz (26). With regard to relation ([33]), this leads to the 



(34) 



following prediction for the solution supported by the single 5-function in the second finite bandgap: 

5 



A 2 V 24 



-£ 2 — V 



As follows from Eqs. (34), inside the second finite bandgap (32) this solution exists for a = —1, i.e., solely for the 
attractive nonlinearity. 
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FIG. 3. The same as in Fig. [2] but for the attractive nonlinearity. In panels (a)-(c), e = —0.5, and in (d)- 



20 



(f), e 



In the first approximation, the norm of the soliton is 

2 (5/24) -v 



N 9 



e 2 (5/576) + (P/6) - v 



,72 



((5/24) -vf 



(5/576) + (z>/6) - v 



r/2 



(35) 



with v defined as per Eq. (32l. It follows from the plot of the norm versus /i, which is shown in Fig.[5]for e = ±1, that 
the corresponding GS family satisfies the VK criterion, dN/d/i < 0. Nevertheless, as well as in the case of the solitons 
in first finite bandgap, which was considered above, the GSs in the second bandgap, supported by the attractive 
nonlinearity, turn out to be unstable at all values of e. 

In the general case, when the (5-function is placed asymmetrically with regard to the OL potential, GSs may exist 
in the second bandgap in the case of the repulsive nonlinearity, and may be stable in that case. A detailed analysis 
of this case is presented in the next section. 

Examples of numerically found profiles of the solitons in the second finite bandgap, together with their analytically 
predicted counterparts [see Eqs. (33) and (34)], are displayed in Fig. [6] The results shown in both Figs. [6] and [5] 
demonstrate that the perturbative approximation is more accurate for smaller and positive values of e, when the 
attractive (5-function is set at a local minimum of the OL potential. 



C. The semi-infinite gap 



For the sake of completeness of the analysis, we also present results of the application of the perturbation theory 
(for small |e|) to solitons supported by the attractive (5-function (a — —1) in the semi-infinite gap, i.e., at fi < 0. In 
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(a) (b) (c) (d) 

FIG. 4. The norm of the solitons in the first finite bandgap, as a function of /i. Examples for the repulsive nonlinearity are 
given for e = 1 (a) and e = 0.5 (b), and the attractive model is demonstrated for e = —1 (c) and e = —0.5 (d). Numerical 
results and the analytical approximation are depicted by solid and dashed lines, respectively. 




FIG. 5. The norm of the analytically predicted (dashed curves) and numerically found (solid curves) solitons in the second 
bandgap as a function of fj,, for e = 1 (a) and e = — 1 (b), in the case of the attractive nonlinearity. 



the zeroth approximation (e = 0), the soliton is given by solution Q. Then, it is easy to find the first-order correction 
to it in the following form: 



1-2(1 



(j>o(x) 



1 



cos (2x) 



|sm(2|,1) 



(36) 



where 4>o(x) is expression Q [to the first order in e, solution 4>o(x) + 4>\{x) satisfies both the linear part of Eq. ( 11 ) 
at x 7^ 0, and the jump condition |19| at x = 0]. To the same order, the norm of this solution is 



N = 1 



£ (5 + 2/j.) 
3(1 -2/z) 2 



(37) 



which demonstrates that the weak OL potential lifts the degeneracy of the soliton family Q. Note that expression 
(37) satisfies the VK stability criterion, 



dN 
dfi 



2e(ll + 2/x) 
3(1 - 2/z) 3 



< 0, 



(38) 



for jU > —5.5 in the case of e > 0, which corresponds to the attractive (5-function set at a local minimum of the OL 
potential, see Eq. (11); the change of the sign of dN/dfj, at /i < —5.5 does not really matter, as dependence (37) is 
virtually flat in that region. 

Numerical results demonstrate that the entire soliton family is stable in the case of e > 0, obeying the VK criterion 
everywhere in the semi-infinite gap (for a further discussion, see the next section) . Figure [7] presents the respective 
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FIG. 6. Typical examples of broad profiles of the solitons in the second finite bandgap, for e — 1 (a)-(c) and e = —1 (d)-(f). 
The respective values of fi are indicated in each panel. Sol id curves depict numerically found profiles, whereas the dashed lines 
represent the analytical prediction given by Eqs. (33 1- (34 1. 



comparison between the analytical approximation and the numerical results. As expected, the prediction is more 
accurate for smaller e, closely approximating even complex soliton profiles near the edge of the gap, see panel (e) in 
Fig. [7] On the other hand, the solitons centered at the local maximum of the OL potential at e < are definitely 
unstable (not shown here in detail). 



IV. NUMERICAL RESULTS FOR THE MODEL WITH THE SINGLE 5-FUNCTION 

In the case when the local nonlinearity is represented by the single S- function in Eq. ([7]), we have found soliton 
modes in the semi-infinite and two lowest finite gaps, and their stability was investigated by means of numerical 
methods. The stationary solutions were constructed applying the Newton-Raphson method to the respective nonlinear 
boundary- value problem. The stability was then examined by considering a perturbed solution to Eq. ([7]), in the form 
of 



4>{x,t) = 4> s {x) 



ge- m + f 



iX't 



(39) 



where 4> s (x) is the stationary solution, functions g and / are eigenmodes of the infinitesimal perturbation, and A is the 
corresponding eigenvalue, which may be complex in the general case. When substituting expression (391 into Eq. ^ 
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FIG. 7. Panels (a) and (d) display the comparison between the numerical and analytically estimated norms (the solid and dashed 
curves, respectively) for the solitons in the semi-infinite gap, for e = 0.1 and 0.5, respectively. Examples of the soliton profiles, 
taken close to the edge of the gap and deeper inside, are displayed in panels (b),(c) and (e),(f), for the cases corresponding 
to (a) and (d), respectively. In (b) and (c), the analytically predicted profiles are very close to their numerical counterparts, 
making the dashed and solid lines virtually indistinguishable. 



and linearizing, one arrives at the following eigenvalue problem: 



aSix-MM*))' -L A/rH/J' (40) 

with L = —fi — (l/2)c£ 2 /dx 2 — ecos(2x) + 2ad(x — £) (<j> s (x)) . This problem can be easily solved using a simple 
finite-difference scheme, the solution being stable if all the eigenvalues are real. 

The so predicted stability or instability was verified by means of direct simulations of the evolution of initially 
perturbed modes. For this purpose, the standard pseudospectral split-step method and the Crank-Nicolson finite- 
difference algorithm were used. In the course of the analysis, the ^-function was set at different positions within half 
a period of the OL, < £ < 7r/2. In the numerical calculations based on the discretization with stepsize Ax, the 
discrete counterpart of the ^-function was defined so that it took a nonzero value, 5 = 1/(2 Ax), at the single point, 
x = t 



A. Solitons in the semi-infinite gap 



Numerical analysis of the solutions in the semi-infinite gap reveals a single soliton family, which exists only in the 
case of the attractive nonlinearity (<r = — 1). A natural result of the stability analysis is that these solitons are stable 
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FIG. 8. (Color online) The norm versus the position of the attractive 5-function, for the soliton family in the semi-infinite 
gap, at e = 5 and fj, = —4. Here and in similar figures shown below, the dashed line depicts the underlying periodic potential 
(rescaled and shifted upward for the clarity of the picture), while stable and unstable solitons correspond to continuous and 
dotted lines, respectively. 



if the attractive 5-function is located at or near the minimum of the potential. An example for e — 5 and /i = —4 is 
displayed in Fig. |HJ where shift £ of the 5-function from the potential minimum, x = 0, takes values within a half of 
the spatial period of the potential (recall the period is n). In the semi- infinite gap, the soliton's stability complies with 
the VK criterion. In particular, the soliton family is completely stable for £ = 0. On the other hand, the family is 
(quite naturally) completely unstable for x = n/2, when the 5-function is set at the point of the potential maximum. 



B. Solitons in the first finite bandgap 



The numerical analysis of solutions in the first bandgap demonstrates that GSs cannot exist simultaneously for the 
attractive and repulsive nonlinearities. If the 5-function is set close to a maximum of the potential (£ = tt/2, for 
e > 0) , it can support a soliton only with the attractive sign of the nonlinearity. Shifting the position of the attractive 
(5-function towards an adjacent minimum of the potential, the soliton ceases to exist at a critical (threshold) value 
of the coord inate, £ = £thr- At the same point, a new soliton appears in the repulsive case, and exists at £ < £thr- 
Figure |9(a)| describes the norm of the GS versus the 5-function's position within a half of the OL period, for e = 5 
and fj, = — 1, which is close to the middle of the first finite bandgap. These results agree with the prediction of the 
perturbation theory obtained for x = 0, according to which the soliton exists only for a = sgn(e), see Eq. (21 ) (recall 
that s > and £ = tt/2 are equivalent to e < and £ = 0). Close to the existence threshold, the amplitude of the 
soliton diverges, while its width remains approximately constant. Representative examples of such GSs are displayed in 

10(c)f(f) demonstrate a stable soliton (supported by the repulsive nonlinearity) featuring 
e, achieved at x = 0.282. The threshold value, £thn varies with e and /i, as shown in 



10 In particular, Figs. 



Fig._ 

an especially high amplitude 

Fig. 11 Specifically, for small values of e and/or small values of /i (close to the lower edge of the first bandgap), the 
soliton-existence region expands in the case of the repulsive nonlinearity. 

The stability analysis demonstrates that, in the case corresponding to Fig. 



9(a) the GSs in the first bandgap, 



supported by the repulsive 5-function, are always stable. On the other hand, in the case of the attractive nonlinearity, 
local unstable eigenmodes exist, with large (Im { A} > 1) purely imaginary eigenvalues, making all the solitons unstable. 
This finding is not surprising because, as stressed above, in the cases of the repulsion and attraction the GSs tend to 
be located, respectively, close to a local minimum or maximum of the periodic potential. 

The stability analysis was also carried out for other values of e and [i. The solitons are always stable under 
the repulsive nonlinearity, while the stability may change in the case of the attraction. Examples are presented in 
Fig. |12(b)f(e) for e = 5 and £ = 0.15 (repulsion) and £ = 0.4, 0.7, 1.1 (attraction). In the latter case, it is seen 
that, for n close to the upper edge of the first bandgap, there are small stability areas for the solitons (conspicuous 



14 



50 
45 
40 
35 
30 
Z 25 
20 
15 
10 
5 














H=-l . 


(b) 















CO 








-—"Repulsive 








.^-—Attractive 


- 

i 










■;, (d) 









0.5 



1.5 



50 
45 
40 
35 
30 
Z 25 
20 
15 
10 
5 



Attractive 



ii=0.7 



Stable 




(b) 

o- 



(O 



(a) 



0.5 1 

\ 

(b) 



1.5 



FIG. 9. (Color online) The norm versus the position of the nonlinearity-modulating 5-function, for solitons in the first finite 
bandgap, with e = 5 and y, = — 1 (a) or n — 0.7 (b), and both signs of the nonlinearity. The dashed-dotted vertical lines 
represent the respective thresholds, £ t hr = 0.287 (a) and £thr = 0.035 (b). As before, continuous lines refer to stable solitons, 
while dotted and dashed-dotted ones correspond to the strong localized and weak oscillating instabilities, respectively. Shapes 
of the solitons corresponding to marked points in (a) are displayed in Fig. 10 The evolution of the solitons marked in (b) is 
displayed in Fig. 13 In panel (b), only the soliton branch corresponding to the attractive nonlinearity appears, as the region 



of £ for the repulsive sign is too small, while the corresponding solitons have a very large norm. 



for £ = 0.4, and extremely small for larger £). For fixed /i = 0.7 and varying position £ of the attractive <5-function, 
a small stability region is found too, as shown in Fig. |9(b)| Although it is small, the existence of the stability area 
for the GSs in the case of the attraction is a remarkable fact, as it is often assumed that all solitons are unstable in 
finite bandgaps, if the nonlinearity is attractive (see, however, Ref. [TS], where stable GSs were found in the case of 
attraction). Apart from the common "localized" instability discussed above (indicated in the stability diagrams by 
dotted lines), there is additional weaker instability which is often referred to as "oscillatory" one (see Refs. [17l[T8] for 
more details). This instability is characterized by complex eigenvalues, A, and delocalized oscillatory eigenfunctions. 
Techniques which can be used for the detection of such instabilities are described in Ref. |18j . Regions corresponding 
to unstable solitons of this type can be clearly seen in Fig. |9(b)| distinguished in this and other figures by a dashed- 



dotted marking. The oscillatory instability can also be found in the case presented in Fig. 12(e) within a tiny area 
between the non-oscillating unstable solitons and the region of stable solitons, which is very small by itself. Figure fl3] 
illustrates the stability and the simulated development of the solitons' instabilities, for the representative cases marked 
in Fig. [9(b)] 

It is relevant to stress that the stability investigation was carried out repeatedly, varying the size of the spatial 
domain and the number of grid points. By doing so, we have checked that the stable solitons exist indeed, not being 
simply a case of a weak instability. 

Looking at Fig. [12] it is easy to conclude that, unlike the semi-infinite gap, the stability of the GSs in the first 
finite bandgap does not obey the VK criterion (cf. Ref. [55] . where the same conclusion was made for models with 
combined linear and nonlinear periodic potentials). Another noteworthy fact is that, when the <5-function is placed 
anywhere except maxima or minima of the potential (x = 0,tt/2), the GS family does not completely fill the first 
bandgap. For instance, for £ = 0.7, the lower cutoff (existence border) for the soliton in the model with attraction is 
/i thr —2.5, while the lower border of the bandgap is at fi = —2.894. 



C. Soliton solutions in the second finite bandgap 



The numerical investigation was also performed for GSs in the second finite bandgap. For the repulsive nonlinearity, 
a single branch of solitons exists in one half of the potential period in the (N, £) plane. On the other hand, in the 
case of the attractive nonlinearity, there are two different branches, one located around the minimum of the potential 
and the other one - around its maximum. A typical example is shown in Fig. 14(a) for e = 5 and /i = 2, close to 
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FIG. 10. Shapes of the solitons in the first finite bandgap, for the repulsive (a)-(b) and attractive (c)-(d) single-5-function 
nonlinearities. The solitons correspond to marked points in Fig. 9(a) with e = 5, fj, = —1 and (a) £ = 0, (b) f = 0.18, (c) 
£ = 0.4 and (d) £ = n/2. The vertical arrow in each panel denotes the position and sign of the 5-function. A noteworthy 
example is the soliton in panel (e), with a particularly large amplitude, obtained at £ = 0.282, which is very close to £thi- Panel 
(f) illustrates the stability of this soliton in direct simulations. 
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FIG. 11. (a) The existence threshold, £thr, for the solitons in the first finite bandgap, as a function of e, for /i = —1. (b) The 
same, but as a function of fj, for fixed e = 5. 
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FIG. 12. (Color online) The norm of the gap solitons versus fi, in the first finite bandgap, for e — 5, with the repulsive 
5-function placed at £ = (a) or £ = 0.15 (b), and the attractive 5-function placed at £ = 0.4 (c), £ = 0.7 (d), £ = 1.1 
(e), and £ = n/2 (f). As above, the stable and unstable solitons correspond to continuous and dotted lines, respectively. A 
soliton subject to an oscillating instability can be found in the case (e), in an extremely narrow region (barely visible on the 
scale of this figure) between the stable and locally unstable sections. Blue vertical stripes on both sides indicate Bloch bands 
between which the first bandgap is sandwiched. The dashed-dotted vertical lines represent thresholds at which the solitons' 
norm diverges: fi t hr = —0.141, —1.572, —2.5, —2.858 for (b), (c), (d) and (e), respectively. 




FIG. 13. The evolution of solitons supported by the single attractive 5-function, for the case shown in Fig. 9(b) Panel (a) 
displays an example of a stable soliton, obtained for £ = 0.4. Unstable solitons, subject to the oscillatory instability or the 
strong localized instability, are exhibited for £ = 0.9 (b) and £ = 1.3 (c), respectively. 
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FIG. 14. (Color online) The same as in Fig. |9j but in the second finite bandgap, for fi = 2 (a) and fj, = 1.4 (b). Note that, 
in very narrow regions of £, such as the left margin in (b), the norm of the (unstable) solitons, supported by the attractive 
instability, is especially large, therefore it cannot be displayed on the scale of this figure. The two thresholds (as explained in 
the text) are represented, as usual, by dashed-dotted vertical lines. 



the middle of the second bandgap. In this bandgap, there are two thresholds in the region of < £ < tt/2. Similar 
to what was seen in the first bandgap, the soliton's amplitude diverges at the threshold, while the soliton's width 
remains finite. 

With the attractive 5-function, GSs are always unstable in the second bandgap (featuring the localized instability) , 
on the contrary to the situation in the first bandgap, where a small stability region was found for the case of attraction, 
see Figs. |9]and[T2] On the other hand, in the case of the repulsive 5-function, the stability changes with the variation of 
parameters £, e and \x. In particular, in the situation displayed in Fig. 14(a) both localized and oscillatory instabilities 
can be found, and two stability regions are present. If the repulsive <5-function is set at £ = 0.5, and e = 5, in which 
case the second bandgap is 1.05 < \x < 3.724, stable solitons are found only at fi close to the lower edge of the gap, 
see Fig. [15] For instance, at — 1.4 the GSs arc stable at almost all values of £, except for a narrow interval where 
the localized instability occurs, as seen in Fig. |14(b)| Similar results were found at other values of e. 

Similar to the situation in the first bandgap, the GS families cover the entire second bandgap only for £ = and 
7r/2. For different values of £, there are regions of /i in which no solitons are present, see Fig. 15 



NUMERICAL RESULTS FOR THE MODEL WITH THE TWO SYMMETRIC ^-FUNCTIONS 



Equation ^ presents a natural extension of the model, which includes two (5-functions symmetrically positioned 
around the potential maxima or minima. We here consider the existence and stability conditions for solitons in this 
model, in the semi-infinite and the first finite gaps. In each case, two settings were explored, with either a potential 
maximum or minimum located exactly at the midpoint between the two (5-functions. The stability analysis was 
carried out using the method outlined in Sec. 



IV with the nonlinearity coefficient S(x — £) in Eq.(40) replaced by 



A. Solitons in the semi-infinite gap 



As in the models with the single (5-function, solitons in the semi-infinite gap exist only for the attractive nonlinearity. 
First, we examined the changes that the solitons undergo with the increase of distance £ of each <5-function from the 
potential minimum located between them, which corresponds to e > in Eq. (JsJ) . 

For small values of £, a region of stable symmetric solitons always exists. Increasing £, we reach a bifurcation point, 
after which the symmetric solutions lose their stability (against non-oscillatory perturbations) and a new asymmetric 
branch emerges, that may be partially stable. An example is shown in Fig. 16(a) for e — 5 and /i = —4. In this case, 
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FIG. 15. (Color online) The norm versus (i, in the second finite bandgap for e = 5, with the repulsive 5-function set at £ = 0.5. 
The norm diverges at the threshold /j, = 3.307, marked by the vertical dashed-dotted curve. 



there is a tiny region of stable asymmetric solitons, abutting on the bifurcation point, which is too small to be visible 
in the figure. Closer to the edge of the semi-infinite gap, the bifurcation occurs at higher values of £. For example, 
f bif ( £ = 2, fi = -1) = 0.777, £ bif (e = 5, = -3) = 0.935, f bif (e = 6,/x = -4) = 0.683, and £ bif (e = 8,/x = -6) = 0.491, 
and in all these cases, the asymmetric modes are completed unstable. On the other hand, the bifurcation occurs at 
smaller values of £ for sets of taken deeper inside gap, i.e., when e is smaller and/or /i is more negative, the 

bifurcation giving rise to a conspicuous stability region for asymmetric states in such cases. Typical examples are 
presented in Fig. 18(a)|[(b) for (e,[i) — (2, —4) and (5, —6). Going still deeper into the semi-infinite gap, the value 



of £ at the bifurcation point keeps decreasing. The corresponding domain of stable asymmetric solutions may not 



necessarily emerge exactly at the bifurcation point, but slightly later. An example is displayed in Fig. 18(c) for 
( E ,/i) = (2,-6). 

Bifurcation diagrams in the ([/,, 9) plane, displayed in Fig. |18(d)] for £ = 0.4 and e — 2 and 5, demonstrate the SSB 
(spontaneous symmetry breaking) of supercritical type. For instance, at e = 2, stable branches of the asymmetric 
modes emerge at the bifurcation point, where the symmetric branch loses its stability. On the other hand, for e — 5 
the asymmetric branches are not immediately stable after the bifurcation point, as the VK criterion is satisfied for 
them only when further decreasing fi. 

At large values of £, the asymmetric soliton gradually transforms into a fundamental one, pinned to either of the two 
(5-functions. As concerns the symmetric modes, with the increase of the distance between the (5-functions (2£), they 
transform into two-soliton bound states which never regain the stability they had prior to the bifurcation. Specifically, 
near the OL minimum points, where the fundamental soliton is stable (see Fig. 16(a) ), eigenvalues accounting for the 
instability of the bound state decrease as its two constituents are pulled farther apart with the increase of £ (but still, 
never fall exactly to zero). 

As in the two-5-functions model without the lattice (e = 0) [51] , antisymmetric states exist too and are unstable 
at small values of £. As £ increases, they develop into antisymmetric two-soliton bound states, for which unstable 
eigenmodes could not be found (featuring instead zero eigenvalues) exactly in the region where the stable fundamental 
soliton is supported by the single ^-function. Direct simulations have shown that, while such bound states are stable 
against symmetric disturbances, asymmetric perturbations break them into mutually incoherent fundamental solitons. 



Figure [17 



Fig. 16(a) 



presents several examples of the solitons of all the aforementioned types, in the case corresponding to 



Similar analysis was carried out for two attractive (5-functions placed on both sides of a local maximum of the 
periodic potential, which implies e < in Eq. ([8]). Figures |l6(b)| and 19 present a typical example of the results in 
the (£, iV) plane, for e = —5 and fi = —4. In this case, both the symmetric and asymmetric solitons are unstable at 
small values of £. In other aspects, the results resemble those reported above for the case of two (5-functions placed 
symmetrically around a potential minimum. 

It is relevant to compare these results with those reported in Ref. [5T] for two attractive (5-functions in the absence 
of the periodic potential (e = 0). In that case, exact analytical solutions are available for the pinned states of 
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FIG. 16. Branches of symmetric and asymmetric soliton modes in the semi-infinite gap of the model with two 5-functions, for 
fi — —4 and e = 5 (a) or —5 (b), which correspond, respectively, to the local minimum or maximum of the periodic potential 
located between the ^-functions. As before, the continuous and dotted lines represent stable and unstable portions of the soliton 



families. Circles correspond to representative examples of solitons that are shown in panels (a)-(d) of Figs. 17 and 19 
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FIG. 17. Typical profiles of modes generated by two attractive (^-functions in the semi-infinite gap. Panels (a)-(d) correspond 
to the points marked by circles in Fig. |16(a) The arrows indicate the position of the 5-functions. Stable symmetric (a) and 
asymmetric (b)-(c) solitons are demonstrated, as well as an unstable symmetric bound state (d). In addition, panels (e) and 
(f) show, respectively, an unstable antisymmetric soliton, and an antisymmetric bound state which is unstable with respect to 
asymmetric perturbations. 
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FIG. 18. (a)-(c) Branches of the symmetric and asymmetric states in the (£,N) plane near the bifurcation point, in the semi- 
infinite gap of the model with two 5-functions, for (a) e = 5, (i = —6, (b) e = 2, fj, = —4, and (c) s = 2, (x = —6. (d) The 
bifurcation diagrams in the (fi, 6) plane, for £ = 0.4 and e = 2 and 5. 
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FIG. 19. Examples of solitons corresponding to the points marked by circles in Fig. 16(b) Panels (a) and (b) demonstrate 
unstable symmetric and asymmetric solitons. An example for a stable asymmetric soliton is displayed in (c). Panel (d) 
shows an unstable symmetric bound state. In addition, examples of an antisymmetric soliton, that features the strong local 
instability, and antisymmetric bound state, which is unstable against asymmetric perturbations, are shown in panels (e) and 
(f), respectively. 



all the types — symmetric, asymmetric, and antisymmetric. The symmetric states are stable before the symmetry- 
breaking bifurcation, and unstable after it, terminating at final TV. The bifurcation is of an "extreme subcritical" 
type, with branches of the asymmetric states going backwards and never turning forward, hence they are completely 
unstable. In fact, these results, obtained with the ideal (5-functions, are degenerate. The numerical analysis with 
regularized (5-functions lifts the degeneracy, demonstrating that the branches of the asymmetric modes eventually turn 
forward, stabilizing themselves. Simultaneously, the family of the symmetric modes (unstable past the bifurcation 
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point) extends to N —> oo. Antisymmetric states are completely unstable in that model (they become stable if the 
regularized 6- functions are made broad enough) . The comparison with the present findings suggests that the addition 
of the periodic potential also lifts the degeneracy of the system, even if the (S-functions are kept in the ideal form. In 
fact, this is similar to how the inclusion of the weak periodic potential lifts the degeneracy of the soliton family Q 
supported by the single attractive (^-function and stabilizes the family, see Eqs. (36) and (37 1. 



B. Solitons in the first finite bandgap 



Similar to the case of the single (^-function, GSs in the first finite bandgap can be found for both the attractive 
and repulsive nonlinearities. We again start by considering the pair of <5-functions placed around a minimum of the 
potential. For the case of two attractive (^-functions, antisymmetric solitons exist in the region of < £ < it — £ t hr, 
where £thr is the same threshold as in the model with the single ^-function. Unlike what was observed in the semi- 
infinite gap, in the first finite bandgap the asymmetric solitons bifurcate from the antisymmetric branch, see an 
example in Fig. 20(c) for e = 5 and £ = 1. In particular, this setting features a closed bifurcation loop, where both the 
direct and the reverse bifurcations are of the supercritical type. In the case presented in Fig. 20(c) there is a small 
region of stable antisymmetric solitons, obtained for high values of fi, near the upper edge of the gap. Decreasing 
fi, the soliton is destabilized by an oscillatory instability, both the antisymmetric and asymmetric branches being 
strongly unstable past the bifurcation point. A typical example, plotted in the (£, iV) plane for e = 5 and ji = — 1, 
20RP) 



is displayed in Figs. 



where a small section of the antisymmetric branch is stable. The stability analysis, 
carried out for different values of e and /i, indicates a somewhat larger stability region for larger e, and for ji taken 
closer to the upper edge of the bandgap. The stability region gradually disappears at smaller values of e, or close to 
the lower edge of the bandgap. A weak oscillatory instability also occurs in the present case. In particular, in the 
situation corresponding to Fig. 20 a very small section (too small to be visible in Fig. |20(a) ) of the corresponding 



weakly unstable antisymmetric solitons exists, starting from the edge of the stable region and extending to slightly 
smaller values of £. 

In this setting, an additional threshold, was found (which is not related to that in the single- J- function model, 
£thn see Fig. 11). The new threshold serves as the upper boundary for always-stable symmetric solitons generated 

Simultaneously, the same threshold is a 
< £ < 7T — £thr, in the case of the attractive 



by a pair of closely placed repulsive (^-functions, as shown in Fig. |20(b) 



ZT2T 



lower border for a branch of unstable symmetric states that exists at £ thr 

nonlinearity. Ty pica l examples of the soliton profiles of the symmetric, antisymmetric and asymmetric types are 



(2) 

displayed in Fig. 21 Further, Fig. 22 shows £\ 1 as a function of e and /i (within the first finite bandgap) 
symmetric GSs are stable in the case of repulsion, and unstable under attraction. 

Modes found at £ > n — £ t hr may be considered as bound states of two fundamental GSs. For both signs of the 
nonlinearity, symmetric and antisymmetric bound states are found precisely where their fundamental counterparts 
This can be seen in Fig. 20(a) comparing the location of the bound-state branches with those obtained in 

cf. 



All the 



exist. This can be seen in Fig. |20(a) 
the model with the single 5-function. 



Fig. 9(a 



It is not surprising that both symmetric and antisymmetric 
bound states are strongly unstable close to the potential maximum, in the case of the attractive nonlinearity. For the 
repulsive nonlinearity, the symmetric bound states experience weak local instability in the regions near the minimum 
of the OL. On the other hand, their antisymmetric counterparts may appear to be stable in terms of the eigenvalues, 
and in direct simulations with respect to symmetric perturbations, but they split into their fundamental constituents 
when asymmetric disturbances are imposed (similar to the antisymmetric bound states in the semi-infinite gap.) 
When the two (^-functions are placed around a maximum of the periodic potential, the picture is somewhat simpler. 

(2) 

In this case, the secondary threshold, £ th] ., does not exist, while, for the attractive nonlinearity, all types of the soliton 
families — symmetric, asymmetric and antisymmetric ones — exist in the region of < £ < n/2 — £ t h r (not shown here). 
Also, on the contrary to the previous setting, the asymmetric solitons bifurcate from the symmetric branch (not from 
the antisymmetric one) , all the solitons being unstable in the case of the attractive nonlinearity. Results for the bound 
states are qualitatively similar to those reported for the minimum-centered configuration, with the difference that, 
in the present case, the weak localized instability occurs for the antisymmetric solutions, while the symmetric bound 
states are the ones which are unstable only against asymmetric perturbations. 



VI. CONCLUSIONS 



In this work, we have introduced two settings that combine the attractive or repulsive nonlinearity, concentrated 
in one or two points, and the linear periodic potential. For the model with the single (5-function, we have found 
stable solutions in the semi-infinite and two lowest finite gaps. In particular, in the case of the attractive nonlinearity 
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FIG. 20. (Color online) (a)-(b) The soliton's norm versus half the distance between the two <5-functions, f , in the first finite 
bandgap, for e = 5 and fi = — 1. In panel (a), the asymmetric branch bifurcates from the antisymmetric one, in the case of the 
attractive nonlinearity. At larger values of £, the soliton branches correspond to unstable antisymmetric bound states of two 
solitons, for either repulsive or attractive nonlinearity. (b) Branches of symmetric states are shown for both the repulsive and 
attractive nonlinearities. Not shown in (b) are the unstable symmetric bound states, for larger £, as this part of the diagram 
is virtually identical to the one in panel (a), for the antisymmetric bound states. The dashed-dotted vertical lines represent 
the existence thresholds, n ± £ t h r = 2.853, 3.43 (periodic with period n) and = 0.588 [only in (b)]. Circles indicate solitons 



whose profiles are displayed in Fig. 21 (c) The bifurcation diagram for the antisymmetric and asymmetric modes, in the (p, I 
plane, at e = 5 and £ = 1. 



(<7 = —1), the degenerate family of exact soliton solutions exists in the absence of the periodic potential, being fully 
unstable. Even a weak potential lifts the degeneracy and stabilizes the entire family in the semi-infinite gap, provided 
that the <5-function is set in a finite region around a local minimum of the periodic potential. The stability of this 
soliton family agrees with the Vakhitov-Kolokolov criterion. 

In the first finite bandgap, GSs (gap solitons) have been found for both attractive and repulsive nonlinearities, 
although they do not coexist: if the S- function is placed in a finite area around a local minimum of the periodic 
potential, the GS exists only under the repulsion, and in the remainder of the period of the potential, GS can be 
supported solely by attractive nonlinearity. In the first bandgap, all the GSs are stable in the case of the repulsion, 
while the soliton pinned to the attractive ^-function is stable only in a small region, if any. 

In the second bandgap, two soliton branches, centered around the attractive (^-function set at either the minimum 
or maximum of the periodic potential, were found. While none of them is stable, stability regions were produced for 
the soliton branch supported by the repulsive (Sa-function. It exists exactly for those values of shift £ of the ^-function 
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Figs 



20(a) (b) 



Examples ol solitons generated by two <5-functions in the first finite bandgap, which correspond to the marked points in 
As in Fig. 10 the arrows indicate the position and sign of the (^-functions. Stable and unstable antisymmetric 
as well as an unstable asymmetric one, are shown in (a), (b) and (c), respectively, in the case of two attractive 
Symmetric solitons, stable and unstable, are shown in panels (d),(e) and (f), for the repulsive and attractive 



solitons 
(^-functions 
nonlinearity, respectively 
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relative to the underlying potential at which both aforementioned branches are absent in the case of the attraction. 

In addition to the numerical results, analytical predictions, based on the perturbation theory for the Mathieu 
equation, were presented for the case when the 5-function is positioned exactly at the maximum or minimum of the 
periodic potential. This approximation was developed for the solitons in the semi-infinite and two lowest finite gaps, 
showing a good accuracy in comparison with the numerical findings, in all the cases. 

In the model with two separated local nonlinearities, represented by the symmetric pair of the 5-functions, the 
numerical analysis was carried out for the semi-infinite and the first finite gaps. Families of symmetric, antisymmetric 
and asymmetric solitons were found, with respect to the symmetric set of the two 5-functions. In the semi-infinite gap, 
with the 5-functions set symmetrically around a minimum of the potential, the symmetric solitons are stable up to the 
symmetry-breaking bifurcation point, from which asymmetric branches emerge, that turn out to be partially stable. 
Antisymmetric modes were found too in the semi-infinite gap, being completely unstable. With the two 5-functions 
placed symmetrically around a point of the potential maximum, all the soliton families obtained at small separations 
between the 5-functions, as well as all the two-soliton bound states, are unstable. 

In the first finite bandgap, asymmetric solitons bifurcate from the antisymmetric branch, if the potential minimum 
is set between the attractive (5-functions. This antisymmetric branch has a very short stability segment, while the 
asymmetric GSs are completely unstable. Under the repulsion, there is a stable branch of symmetric modes centered 
around the minimum of the potential. At a certain threshold value of the separation 2£ between the 5-functions (which 
is different from the threshold found in the single-5-function setting) , the symmetric mode switches from the repulsive 
to attractive nonlinearity, simultaneously loosing its stability. Finally, in the configuration with the local maximum 
of the periodic potential fixed at the midpoint, no stable solitons were found for sufficiently small separations 2£. For 
both configurations, with the midpoint coinciding with either the maximum or minimum of the underlying potential, 
all the families of two-soliton bound states are unstable against various perturbation modes. 

This work may be naturally extended in other directions. In particular, it may be interesting to consider the 
nonlinearity represented by a periodic array of 5-functions. A challenging issue is to analyze similar models in two 
dimensions. 
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